IFJPAN-V-04-09 



Hierarchically Organized Iterative Solutions of the 
Evolution Equations in QCD* 



S. Jadach, M. Skrzypek and Z. Was 

Institute of Nuclear Physics, Polish Academy of Sciences, 
ul. Radzikowskiego 152, 31-342 Cracow, Poland, 



Abstract 

The task of Monte Carlo simulation of the evolution of the parton distributions in 
QCD and of constructing new parton shower Monte Carlo algorithms requires new 
way of organizing solutions of the QCD evolution equations, in which quark— gluon 
transitions on one hand and quark— quark or gluon— gluon transitions (pure gluon- 
strahlung) on the other hand, are treated separately and differently. This requires 
certain reorganization of the iterative solutions of the QCD evolution equations and 
leads to what we refer to as a hierarchic iterative solutions of the evolution equa- 
tions. We present three formal derivations of such a solution. Results presented 
here are already used in the other recent works to formulate new MC algorithms 
for the parton-shower-like implementations of the QCD evolution equations. They 
are primarily of the non-Markovian type. However, such a solution can be used for 
the Mar kovian- type MCs as well. We also comment briefly on the relation of the 
presented formalism to similar methods used in other branches of physics. 
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1 Introduction to the problem 



The standard QCD evolution equation 
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describes response of the parton distribution function (PDF) D k (t,x) to a change of 
the large energy scale Q = exp(i). Variable x is identified as a fraction of the hadron 
momentum carried by parton of the type k (quark, gluon). Evolution kernel 7 is calculable 
within perturbative QCD. 

The above evolution equation (CQ) is an important ingredient in many QCD pertur- 
bative calculations. It can be solved using variety of the numerical methods, including 
Monte Carlo method. The knowledge of D k (to,x) at certain initial to, is required for 
solving evolution equation at other t > to- The initial PDF is fitted to experimental data. 

The principal aim of this work is to derive the followinganalytical solution of the above 
QCD evolution equations 
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where we have isolated flavour conserving (bremsstrahlung) part of the kernefl = 
Skfikk from the total kernel, J'fcj = T^- + IPjy, and G^ fc (ti,t ,z) is the solution of the 
following simplified evolution equation, similar to eq. (pO), 

dtG% k (t, t , z) = yf fc (t, •) (g) Gf fc (t, t , •) (z), 



(3) 
See 



see Section E3] for the details. The boundary condition is G kk (to,t , z) = 5(1 — ; 
also fig. Q] for graphical representation. 

It is important to provide formal proof of eq. (J2J), because it is a critical ingredient 
in several new Monte Carlo algorithms of the non-Markovian type3 described in refs. [1] 
and [2], and possibly in other future works. 



^ince T^ fc = 0, only flavour changing indices fej ^ fcj_i enter in the flavour sum in eq. ©• 



It can serve as basis of a novel type of the Markovian MC as well. 
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Figure 1: Kinematics in the 2-level emission tree of eq. (|2])- 

Let us now explain in details notation used in eqs. ([H13]). In function D k (t,x) variable 
1 > x > is the fraction of the hadron momentum carried by the parton of the type 
k = G,qi,qi, i.e. gluon, quark or antiquark, at the high energy scale Q, conveniently 
translated into the "evolution time" variable t — In Q. In QCD the PDF represents the 
wave function of the hadron close to the light-cone. See ref. [3] for an expert discussion 
on the precise meaning of PDF in QCD, in a wide context of the so-called factorization 
theorems [U, [5J [6] in the gauge Quantum Field Theories. 

In this work we shall restrict ourselves to the most common QCD evolution equations of 
the DGLAP type [7j, with the kernel splitting functional incorporating the QCD coupling 
constant (for the sake of the simplicity of notation) 

9 kj (t,z) = ^-P kj {t,z). (4) 

7T 

The QCD kernel functions are singular, with singularities of the type (ln(l — z) n / (l — z)) + . 
We shall typically regularize them with the help of an explicit small infrared (IR) cutoff 
parametecl e as follows: 

? kj (t, z) = --J> s kk (t, e)S kj 5(l -z) + 9%(t, z), 9%(t, z) = ? kj (t, z)Q{l -z-e). (5) 

The important Sudakov formfactor to) is directly related to the virtual part of the 
kernels: 

i 

$,(Mo) = J dt'? 5 kk (t',e). (6) 

to 



3 The DGLAP kernels in the MS scheme were calculated in QCD at the two levels beyond the leading- 
logarithmic (LL) approximation. 

The infinitesimal parameter e can be ^-dependent, without any loss of generality in the following 
treatment. 
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Finally the (bremsstrahlung-type) auxiliary distribution 

G* k (t,t ,z) = 5(l-z) 

t i 

r n 
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is an iterative^] solution of the flavour-diagonal evolution equation of eq. ([3]). See also 
fig. [2] for graphical representation of the above gluonstrahlung segment of the evolution. 
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Figure 2: Kinematics in the pure bremsstrahlung emission tree of eq. (pp). 



2 The solutions 

If our only aim was to prove the correctness of eq. (j2j) as a solution of eq. ([T]), then 
the simplest approach would be just to substitute it into this equation and check with a 
little bit of algebra that indeed it is the solution. Our aims are however more general: 
(i) to derive eq. ([2]) in a more systematic way, (ii) to understand better its relation to 
the other widely known and used iterative solutions of eq. (CEJ), (iii) to prove that its 
exclusive content, in terms of the fully differential distribution in all variables and Zi, 
i — 1,2, ...,n, for each n, is exactly the same as in other iterative solutions, commonly 
used in the MC approaches. 

Having all the above in mind, let us proceed methodically, first with deriving solution 
of the evolution of eq. ([T]), in terms of a time-ordered exponential, widely used in the 
literature. Next, we shall present first example of the derivation of eq. §Z§ by means of re- 
organizing the evolution equation and solving it once again. Then, we shall present second 
example of the derivation, in which the above time-ordered exponential is algebraically 
reorganized (transformed) into eq. ([2]). Finally the third derivation of eq. (T5]) based on 
straightforward reorganization of the multiple sums and integrals will be included in the 
Appendix. 



3 We shall explain in the next Section why we call it "iterative" . 
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2.1 Time-ordered exponential 

The solution of eq. ([I]) can be established quickly and rigorously, for instance by means 
iteration, as a time-ordered exponential of the kernel operator D 5 in the vector (linear) 
space indexed by one continuous variable x and one discrete k. More precisely, eq. §T§ in 
a more compact matrix notation reads 

d t D(t) = P(t) D(t) (8) 

and its solution in the same compact matrix notation is given by 

D(t) = exp Q" P(*')^ D(t ) = Gp(t,t )D(t ), (9) 

where we employ the following well known time- ordered exponential evolution operate^ 

G H (t,to) = G(H;t,to)=exp( \ B.(t')dt , \ = 1 + ^11/ dt A>t^MU), (10) 



to 



n=l i=l 



to 



for t > to- It is familiar to all readers from the textbooks of the Quantum Mechanic^. 
We define function 9 x>y to be equal 1 when x > y and equal otherwise. 

The compact solution of equation (Q, can be translated into more traditional integro- 
tensorial notation, with explicit sums and integrals: 
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where 8 



x=y 



5(x — y). We realize, that the above formula, although rigorous and very 
elegant, is useless for any practical evaluations, because it contains polynomials of the 
negative and singular terms comming from products of — 7 s 5(1 — z) factors. In principle 
these inconvenient terms can be resummed (exponentiated) with the help of the direct but 
tedious algebra on the multiple flavour indices and z-integrala^]. In the following section 
we shall do it using more elegant methods. 



6 Here and in the following we define Yl7=i = An-An-i ■ ■ ■ ^2^1 ■ 

7 However, we do not require H to be hermitian and G to be unitary. 

8 Similarly as the method used in the Appendix to re-sum another part of the kernel. 
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Having denned the time ordered exponential evolution operator G#, let us quote its 
basic features and extend its definition for the latter use. The well known ru 

G H (t, t x )G H (t x , t) = G H (t, t ), t>t x > t (12) 

helps to manipulate products of the time-ordered exponents. We may also define the 
inverse operator for the "backward evolution" (t < to) as follows 

G H (t ,t) = G^(t,t ), t<t , (13) 

where the inverse operator is constructed usin Gtf(t,t ) = G ( _H)(t,t ). With help of 
the above definition validity of eq. ([121) can be extended to any t x . 

2.2 Derivation by reorganizing evolution equation 

In the following we show how to resume singular 7 s terms by going back to the evolution 
equation, reorganizing it and solving it once again. We are going to show this standard 
trick in a detail, because, subsequently, we shall generalize it to the case of an arbitrary 
part of the kernel (instead of the T" 5 part). It is essentially a warm-up example. 

2.3 Resuming virtual part of the kernel — warm-up example 

Inserting explicitly regularized kernel, our evolution equation takes the following form 

d t D k (t,x) = -9{ k {t) D k (t,x) + J2y%^-)®D 3 (t,-)(x). (14) 

j 

It can be transformed int d_j] 

d t (e* k ^D k (t,x)\ = ^e* fc( *'' o) ?®(t,-)e-^ ( * A)) ® e*^> to) D^t, -)(x). (15) 

3 

Changing slightly notation the above is transformed into 

d t D k (t,x) = (t, •) ® D^t, -)(x), 

j 

D k (t, x) = exp($ fc (t, to))D k (t, x), ^ 
■?%(t, z) = exp($ fe (t, t Q )P%(t, z) exp(-$,(t, t )) 

or in an equivalent compact matrix formulation it reads 

d t b{t) = p (t) bit). (17) 

9 Chapman-Kolmogorov-Smoluchowski-Einstein relation, see ref. [H1H]. 

10 The algebraic proof that G^{t, to)Gii(t, <o) = I, using eq. ([TO]) , we leave to the reader. The matrix 
elements of G _1 can be non-positive and highly singular. 

11 This integro-differential form is exposed in the QCD textbooks, see ref. [TU], and is also used in the 
numerical evaluation (evolution) of PDFs using non-Monte-Carlo methods, for instance ref 
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The time ordered solution 



D(t)=exp( I P e (t')dt'\ D(t ) = Gpe(t,t ) D{t 



to 



;is) 



of the evolution equation is widely known and exploited routinely in many practical 
evaluations of solutions of the QCD evolution. It is usually written in the traditional 
integro-tensorial representation similarly as eq. (TIT]) , in terms of of the initial D(t ) and 
the product of CP 8 , taking the following familiar shape: 



OO r- n „ „ 
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dx l[V% k ._ 1 (t i ,z i )e-^-^ ti ' ti -^ D ko (t ,x )6(x-x l[z i ). 
i L i=i -I i=i 



(19) 



Trivial identity 



e -*i(*i.*o) e *i(*i-i,to) 



(20) 

was also employed. The above solution of the evolution equation is used as a basic 
formula in the Monte Carlo evaluation of the PDFs using Markovian MC algorithms, see 
for example ref. [2]. We shall refer to this solution as non-hierarchic iterative solution of 
the evolution equation. 

Let us remind the reader, that in the above warm-up exercise we have resummed 
the relatively simple component 7^ k (t, z) = — 5jfc(P£ fc (£)5(l — z) of the evolution kernel, 
which was completely diagonal, both in k and in z. In this special case G^i is trivially 
calculable, contrary to more general case of non-diagonal y B discussed in the following. 

Let us now come back to our principal aim, proving eq. (j2J), where less trivial compo- 
nent of the kernel will be isolated/resummed. 



2.4 Resumming gluonstrahlung — the real thing 

In order to prove eq. ([21 , we need to resum (exponentiate) the following part of the kernel 

?f k (t, z) = 8 jk ? kk (t, z) = -5 jk ? s kk (t)6(l -z) + 5 jk ?® k (t, z), (21) 

which is diagonal in the flavour indices, but not in z. This part of the kernel is always 
IR divergent and generates multiple gluon emission process, that is gluonstrahlung. The 
remaining flavour- changing part of the full kernel is defined as 7 A = "J , — 'J )B . The original 
full evolution equation and its solution read 

«9D(t) = (p A (t) + P B (t))D(t), D(t)=exp^ (p A (t') + P B (t')yt'^ D(t ). (22) 
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At this point, the G^-function of eq. (JTj) can be identified with the following operator 

Gb(Mo) = G pB {t,t ) = exp Q" P B (t')dt'^J , (23) 

where 



d t G B (t,t ) = P B (t)G B (t,t ), (24) 

see also eq. ([3]). 

In order to derive eq. (T5]) we proceed analogously as in the derivation of eq. (TT9l ; 
we shall introduce GA(t,t ) in the evolution equation, similarly as we have introduced 
exp ($fc(t, to)) • The main complication will be in the non-commutative nature of G^(t, t ). 
Let us introduce in the evolution an equation auxiliary PDF 



D(t) = G^Mo) D(t), D(t) = G B (t,t ) D(t), 
getting after the differentiation 

G A (t,t ) dD(t) + (dG A {t,to))D{t) = (p A (t) + P B (t)) G A {t,t ) D(t). 



After inserting eq. (J24I) we obtain 

G B (t,t ) <9D(t) + P B (t) G B (t,t )f>(t) = (p A (t) + P B (t)) G A (t,t ) D(t). 
The term proportional to P gets eliminated 

Gfl(t,to) «9D(t) = P A (t) Gb(Mo) D(t). 
and we return to the usual evolution equation 

dD(t) = P A (t) D(t), P A (t) = G^Mo) P A (t) G B (t,to) 
with the usual solution 

D(t) = G B (t,t ) exp ^P A (t')dt^ D(t ). 



(25) 



(26) 



(27) 



(28) 



(29) 



(30) 



The last step on the way to eq. (j2J) is elimination of the operator G^ 1 (t, to) being part 
of P. The reason for that is that G _1 is not well suited for any numerical evaluation, 
especially of the MC type, due to alternating sign in the exponential expansion, hence it 
is better to eliminate it from the final result. It is done with the help of the following 
identity 



G B (t,t Q ) exp / P A (t')dt 



to 



G B (t,t Q ) + J2 



n=l 



I! / dtiOt^t^Gsit^t^iU 



(31) 



G B (t\, tc 



where t n+ i = t. This identity is derived rather easily by inspecting each n-th term in the 
expansion of the time ordered exponent and applying the following relational (analogous 
to eq. (JUD) 



Gb(£i+1)£o) G B 1 (tj,to) — Gb(U+i,U) Gs(ti,to) G B 1 (ti,to) — Gb^i+i,^) 
for each pair G^G^ 1 sandwiched between adjacent P^'s. The final solution reads 
D(t) = G B (t,t ) D(t )+ 



(32) 



71=1 



" pt 

I [ / ^>ti_iGs(t«+i, tj)P A (tj 
i=i 



G B (ti,t ) D(t ). 



(33) 



When translated into the integro-tensorial notation, the above formula turns out to be 
identical with our target eq. (TSJ). In this way we have completed its derivation. 

It is now obvious why eq. (J2]) we call a hierarchic solution of the evolution equation. 
It is because its components G^ are solutions of another simpler evolution equation 
(gluonstrahlung) of its own. Higher level evolution embeds lower level simpler evolution 
as a building block. 



2.5 Derivation by reorganizing time-ordered exponential 

A disadvantage of the derivation presented above is that it exploits the inverse evolution 
operator G -1 , which is in a general case difficult to define properly, while it drops out 
from the final result of eqs. ([2]) or (|33|) anyway. The natural question is therefore whether 
we could derive eq. (|2J) without introducing the operator G -1 in the intermediate stages 
of the proof. 

Furthermore, going back to the modified evolution equation and solving it once again 
obscures the relation 13 ! between variables (k%, Zi) in the non-hierarchic solution of eq. ( 1191) 
on one hand and the hierarchic one of eq. (j2J) on the other hand. In the following we shall, 
therefore, present an alternative example of the derivation of eq. (Q without explicit use 
of the inverse evolution operator G _1 . In such a case, the relation between variables 
(ki,Zi) in eq. (1TTJT) and eq. (|2J) can be traced back (recovered) more easily. 

The following derivation will be strongly reminiscent to a derivation of identity F(x, y) = 
exp(x + y) = exp(x) exp(y) by means of the Taylor expansion^ with respect y i.e. 
n^v) = En=o(x n /nl)d2F(x,y)\ x=0 . 

Let us introduce slightly modified evolution operator 

G' H (t,t ) = G'(H;Mq) = G(H;M ) 9 t > to , (34) 

where G was already defined as the time-ordered exponential in eq. fflOl) . The additional 
^-factor ensuring t > to is wm make the following algebra more compact. We define 

12 See eqs. (|12M13|) and the accompanying discussion. 
13 Such a relation is relevant for parton shower applications. 

14 Note that such a derivation is almost equivalent to a direct multiplication of the infinite sums for 
exp(x) and exp(y), but more transparent algebraically. 
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H = Ha = B + A A (we shall set A = 1 at the very end of calculation). The whole 
derivation relies on the following identit\f^l 

(35) 



<9 a G'h(Mo)= / dti G' H (t,ti) A(ti) G' n (ti,t ), 

J t 



which can be derived using definition of eq. ( flOl) . and reorganizing all integrations over 
ti's. The second derivative follows trivially: 

dlG' H (t,t )= [ dh dt 2 G' H (t,t 2 ) A(t 2 ) G'nfah) A(h) Gk(*i,*o) 

J t 

+ f dt x dt 2 G' H (t,ti) A(t 2 ) G' H (t l ,t 2 ) A(tO G' H (t 2 ,t ) (36) 
= 2! I dt 2 dh G' n (t,t 2 ) A(t 2 ) G' H (t 2 ,ti) A(ti) G' H (ti,t ), 

Jto 



'to 

and the n—th derivative i^ 1 ' 1 



d^G' Ux (t,t ) = n\ f[ ^£dU G' Hx (t i+1 ,U) A(U) ^G^Mo), (37) 

where t n = t. Now, let us use Taylor expansion 

G^Mo) = G' Hx Jt,t ) + -r d n x G' Ux (t,t ) . (38) 

o 

Noticing that G' Hx (t i+ i,ti)\ x=Q = G B (t i+1 ,t;), we obtain 

G / H A (t,to) = G' B (t,t ) + ^A" J] / dUG'^ Xt U)A{U)\G^{t u t Q ), 

n=l i=l V *° / 



(39) 



We may set A = 1 at this point. 

Identifying A = P A , G' n = G' pA+pB and G' B = G' pB we obtain more familiar iden- 



titv 171 



G' pA+pB (t,t ) = + (/^ G'flfe+i^) P^fe)) W^to), (40) 

which leads immediately to eqs. (Q and (j2J). In this way we have completed the second 
proof of eq. ([2]) - this time without any reference to backward evolution operator G -1 . 



15 Quitc similar identity holds for an arbitrary A-dependence in H(A). 
16 Strictly speaking we should use mathematical induction over n to verify this. 
17 As previously we define G' x (ti,tj) — Gx(ti,tj)6 ti>tj . 
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2.6 Straightforward derivation 

In addition to two elegant proofs of eq. (j2J) presented in the previous Sections, we include 
in Appendix third proof, which relies on a rather straightforward method - it starts from 
eq. fl33|) and through tedious reorganization of the sums over flavour indices (change of the 
summation order) and relabeling of the variables transforms it into eq. (j2j). The advantage 
of this third proof is that relation between integration and summation variables in both 
formulas is exposed in a manifest way. This might be useful in the construction of the 
exclusive MC model of the parton shower type. 

3 Discussion 

We are fully aware, of course, that all three derivations of eq. ([2]), shown in this work 
represent a well established mathematical formalism, very similar to that in use in the 
Quantum Mechanics, theory of Markovian processes and renormalization group in the 
Quantum Field Theory. We did not add much to the development of the corresponding 
area of mathematics. Rather, our main aim was to customize this known formalism to 
the specific needs of solving the QCD evolution (also numerically), such that solution of 
eq. (j2J) and the other similar ones are obtained in an effortless and rigorous way. Having 
all this in mind, let us comment on certain selected aspects of the presented formalism, 
on their possible refinements, extensions and applications. We shall concentrate mainly 
on two points: 

• Extension to beyond-DGLAP evolutions in QCD, like CCFM and others. 

• Possible application in the Markovian MCs and the related question of the momen- 
tum sum rules and normalization of PDFs. 

3.1 Extensions beyond DGLAP evolution 

In our definitions of the evolution of PDFs eqs. ([TJ{2]) and the rest of the paper we have 
restricted ourself to DGLAP type [7] evolution, leading-logarithmic (LL) version or its 
next-to-LL extensions. This restriction is however inessential and the validity of our 
derivations can be extended to a more general evolution equation 



in which the dependency of the generalized kernel %(t, x, u) is more general than only 
through the ratio z = x/u. The above more general evolution equation is used for instance 
in the CCFM-type models@of PDF [12]. The DGLAP case of eq. (J2J is obviously covered 

18 It is also closer to the spirit of of the parton shower MC and unintegrated PDFs. 




(41) 
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by eq. (HIT) , with the following identification 



/ \ 1 ^ / x\ a(t) 1 / x\ . , 

3C^(t, it) = -Pjy (*,-) = — U, - • 42 

X V It/ 7T X \ 11/ 

The compact matrix notation used in the time-ordered exponentials can easily accommo- 
date multiplications of the kernels %(t, x, u) such that all relevant algebra in the previous 
sections remains unchanged. Let us only indicate how the product of two kernels gets 
redefined 

(P{t 2 )P{t 1 )) k .{x,u) = J2 [ du'X kr (t % ,x,u^X rj (t u u',u). (43) 

j> Jx 

The reader can easily verify that the rest of the compact matrix algebra in our derivations 
remains unchanged. 



3.2 Sum rules and Markovianization 

As already mentioned, results of this work were instrumental in the modelling QCD 
evolution using non-Markovian type Monte Carlo techniques in refs. [13] and [TJ]. The 
corresponding MC programs simulate DGLAP and CCFM class evolutions. 

However, the solution eq. ([2]) may be also used to construct an interesting example of 
the Markovian MC in which single step in the Markovian chain is a Markovian process 
of its own. Without going into fine details, let us indicate how this can be done. To this 
end we have to invoke momentum sum rule 

I dx xD k {t,x) = (44) 
k Jo 

and reorganize slightly eq. (121 . Staying for simplicity with the DGLAP case (LL and 
beyond) the above sum rule determines virtual part of the kernel 



k J ° 



(45) 



The above is used to set up properly Markovian MC and in particular to split Sudakov 
formfactor into bremsstrahlung part and the rest (flavour changing part) 



$f(t,t )= / dt' [ dzz?t k (i',z). 

Jt Jo (46) 

*£(t,ib)= fdt> I dzzV%(t',z). 
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By means of pulling out Q k (i) and multiplying both sides of eq. (T5J) by x, we obtain 
the following formula suitable for a Markovian MC 



/oo „ 
dz dx U kk (t, t , z') x D k (t , x )5(x - z'x ) + ^ dx 

n=l l 





fc„_i,...,fci,fc 



3=1 



to 



1 


r 1 "i 

n „ 


J dz n+l 







. <=i { 



■ft 



i=l 



(n n+1 \ 
£ - %0 Yl Z i II 2 i ) ' k n = k, 
i=l i=l ' 



where 



U k B k (t 1 ,t ,z)=e^ t ^zG* k (t 1 ,t Q ,z) 



(47) 



(48 



obeys evolution equation of eq. ([7]), with the substitution $fc(tj, U-i) ~> $ k (tiiti-i) an d 
with both side multiplied by z. The evolution operator U obeys nice "unitarity" rule 



J dzU k 3 k (t,t ,z) = l 



(49) 



for any t, t > t , in addition to the usual boundary condition U kk (t , t , z) = 5(1 — z). 

Eq. ( 1471) can be now used to define hierarchic (nested) Markovian Monte Carlo al- 
gorithm. The normalized probability distribution of the forward Markovian step in the 
flavour- changing upper level Markovian process reads: 



bj(tii ki \ti— 1, Xi— i, ki—\) — Ul(ti, ZiZ^Xi—i, ki\ti—i, Xi—i, i) 

= (1 - S^e-^^ ZiVj^^Zi) zpt^U, z'lU.,), 

POO rl rl 

/ dU ^\ dZi / dz 'i u {th z i^i x i-lih\ti-\,Xi-l,h-x) = 1- 

ft, 7 . 



(50) 



For the lower level bremsstrahlung process one may use standard Markovian MC technique 
of ref. pQ. 

Let us discuss selected details of the above scenario. Here, all ki and £j, i = 1,2, ...n 
can be generated before generation of any ^-variables in a separate Markovian algorithm 
with the stopping rule being the usual condition t n +i > t. (See ref. [1] for more details 
on the the Markovian class MC algorithms.) Variables Zi of the flavour- changing kernels 
can also be generated at this early stage. 
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The interesting question is: how and when do we generate z[ according to gluon- 
strahlung operator llf?. _ ifcj (U, If we have known an analytical (even approxi- 

mate) representation of this function, or its precise value from the look-up tables, then 
we could readily generate them before entering into MC simulation of the bremsstrahlung 
subprocess. That would lead us to the use of the constrained MC of refs. [131 EE] f° r the 
bremsstrahlung segments. Alternatively, z[ may come out from a separate Markovian MC 
module simulating gluonstrahlung sub-process starting at U-\ and stopping at £$, with the 
normalized probability distribution of single Markovian step defined in ref . [lj . In this lat- 
ter case we would have simulated in the MC a hierarchic system of Markovian processes, 
with the master flavour-changing Markovian process and many Markovian subprocesses, 
each of them implementing pure bremsstrahlung, flavour conserving, emissions. 

The above hierarchic Markovian MC scheme, although quite interesting, seems to have 
no immediate practical importance. However, it may find applications in some future 
works. 

4 Summary 

The basic aim of this paper is to provide solid technical foundation to other works in 
the area of the Monte Carlo simulation of the evolution of PDFs and parton shower 
in QCD. Our basic result is the solution of the evolution equation of eq. (j2J), which 
was proved algebraically using three method. Its primary application is construction of 
the constrained MC algorithms of ref. [13]. In addition, we also describe hypothetical 
application of such a solution in the Markovian MC algorithm. Although we are aware 
of many interesting relation of the discussed problems and solutions to other areas in 
physics, we did not attempt to elaborate on that too much, in order to keep the paper 
compact and transparent. Let us mention also, that our solution can be used many times 
leading to a nested structure with several levels of the hierarchy. 
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Appendix 



A Combinatorial proof 

We are going to show how to transform non-hierarchic solution in eq. (1T91) (with resummed 
virtual corrections) into hierarchic solution in eq. (j2J) (with resummed gluonstrahlung) 
using straightforward method of changing summation order and relabeling integration 
variables. 

The critical point in isolating two levels in the evolution, flavour- changing transitions 
and gluonstrahlung, will be the change of the summation order in eq (j2D, such that 
one is able to resum separately the pure bremsstrahlung segments obeying ki = 
These segments will form (many) functions G kk , as defined in eq. (J7|). The corresponding 
transformation of the summation order (indexing) looks schematically as follows 

oo oo 

^ ] ^ ] tk n k n -i...kiko — ^ ] ^ ] 

n=0 fc n _i...,fcifc n=0 fcn-i-.-.fcifco 

fcn^fc n _l#---^fcl^fcO (51) 



E 



t 



k (jn) A- (2) fr (1) A- (Jn - l) k (2) k (1) k Ul) k {2) k W k (jo) k (2) k (1) 
jn ,jr„-l...J0 = l 

where we have = ■ ■ ■ = k^ = k^ and the purpose of the upper index in this context 
is simply to show that the same index k is repeated j r times. On the other hand variables 
zf 1 ^ and Tr , r — 1, 2, . . . ,n, m — 1, 2, . . . ,j r are truly different (independent), with 
the upper index truly differentiating them. 

The aim is now to show that one can factorize out the functions G kk and identify 
precisely the remaining functions and integrations. Employing the above index transfor- 
mation in the product of the T-functions we obtain 

{7k nkn (t\ t ] ) ■ ■ ■ 3W4 n) , 4 n) )}ow 1 (4 n) , 4 n) ) 

? fclfcl (4 2) ,4 2) )}y fc2fcl (4 2) ,zj 2) ) (52) 

• • •^ lfcl (4 1) ,4 1) )}^ lfc0 (4 1) ,4 1) ) 

{^P 'k k Q {t^jJ , Zj Q ) . . . ^kokoiti ) 4 ) -P&ofco^l ^ 2 1 )}) 

where curly bracket embrace the diagonal elements ^kk, to be collected into the G kk - 
functions; the remaining, nondiagonal ones, are now clearly isolated. 

Each ^Pki,ki-D ki 7^ ki-± in eq. f)19p is accompanied by an exponential factors. All of 
them (including the first one which does not belong to any CP) are now reorganized as 
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follows: 



<J>i (t (1) t W ) ft, (*« f^l 5> fc i (0) ~) ^ 53 '' 

e ^^K..e^o(i 0) A 0) ) e ^ (t{ \t ) 

=e ^ n (tA n) ) e * fewil (4 n) 4 n,1) ) ^(W) e^o(4 1) ,*o) ; 

where all factors entering into products in (^-functions are shown inside the curly brack- 
ets. Together with the flavour-changing 3"s, the above exponential form-factors look as 
follows: 



The other diagonal 5"s will enter into G^ fe -functions. 

In this rather sketchy way we have shown that indeed eq. (1191) can be transformed into 
eq. ([2]) by means of the straightforward reorganization of multiple sums and integrals. 
(More systematic proof would require completing mathematical induction with respect to 
n.) 

References 

[1] K. Golec-Biernat, S. Jadach, W. Placzek, and M. Skrzypek, Acta Phys. Polon. B37 



(2006) 1785-1832, hep-ph/0603031 



[2] S. Jadach and M. Skrzypek, Acta Phys. Polon. B35 (2004) 745-756, 
Pph/0312355| 



[3] J. C. Collins, Acta Phys. Polon. B34 (2003) 3103, hep-ph/0304122 

[4] J. C. Collins, D. E. Soper, and G. Sterman, Nucl. Phys. B250 (1985) 199. 

[5] J. C. Collins and D. E. Soper, Nucl. Phys. B193 (1981) 381. 

[6] G. T. Bodwin, Phys. Rev. D31 (1985) 2616. 

[7] L.N. Lipatov, Sov. J. Nucl. Phys. 20 (1975) 95; 

V.N. Gribov and L.N. Lipatov, Sov. J. Nucl. Phys. 15 (1972) 438; 
G. Altarelli and G. Parisi, Nucl. Phys. 126 (1977) 298; 
Yu. L. Dokshitzer, Sov. Phys. JETP 46 (1977) 64. 



15 



N. G. van Kampen, Stochastic Processes in Physics and Chemistry. North Holland, 
1981. 

[9] A. L. Fetter and J. D. Walecka, Quantum Theory of Many- Particle Systems. McGraw- 
Hill Book Company, 1971. 

[10] R. Ellis, W. Stirling, and B. Webber, QCD and Collider Physics. Cambridge Uni- 
versity Press, 1996. 



[11] M. Botje, ZEUS Note 97-066, http://www.nikhef.nl~h24/qcdcode/ 



[12] M. Ciafaloni, Nucl. Phys. B296 (1988) 49; 

S. Catani, F. Fiorani and G. Marchesini, Phys. Lett. B234 339, Nucl. Phys. B336 
(1990) 18; 

G. Marchesini, Nucl. Phys. B445 (1995) 49. 

[13] S. Jadach and M. Skrzypek, Comput. Phys. Commun. 175 (2006) 511-527, 
|hep-ph/0504263[ 

[14] S. Jadach and M. Skrzypek, Report CERN-PH-TH/2005-146, IFJPAN-V-05- 

09, Contribution to the HERA-LHC Workshop, CERN-DESY, 2004-2005, 



http://www. desy. de/^heralhc/, hep-ph/ 0509178 , 



[15] S. Jadach and M. Skrzypek, Nucl. Phys. Proc. Suppl. 135 (2004) 338-341. 



16 



